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We review the use of the path integral Monte Carlo (PIMC) methodology to the 
study of finite-size quantum clusters, with particular emphasis on recent appli- 
cations to pure and impurity-doped ^He clusters. We describe the principles of 
PIMC, the use of the multilevel Metropolis method for sampling particle permu- 
tations, and the methods used to accurately incorporate anisotropic molecule- 
helium interactions into the path integral scheme. Applications to spectroscopic 
studies of embedded atoms and molecules are summarized, with discussion of the 
new concepts of local and nanoscale superfluidity that have been generated by 
recent PIMC studies of the impurity-doped "^He clusters. 



1. Introduction 

Over the past 15 years, the path integral Monte Carlo (PIMC) method has evolved 
into a uniquely powerful computational tool for the study of bulk and finite quantum 
systems. In PIMC, one is interested in computing the thermal average of a quantum 
observable O at a given temperature T, which can be expressed with respect to the 
thermal density matrix p{R,R';(3) = {R'\e-'^"\R): 

(O) = Z-'^ j dRdR' p{R, i?'; 0){R\6\R'), (1) 

where R = (ri,r2, . . . ,'rj\j) is a point in the SA^— dimensional configuration space 
of an iV-particle system, H is the Hamiltonian, and f3 = l/ksT. Here Z = 
J dR p{R, R; P) is the partition function. The multidimensional integral of Eq. 
can in principle be evaluated by standard Monte Carlo integration schemes, i.e. by 
taking an average of {R\0\R') over the configurations {R,R'} sampled from the 
probability distribution Z~^ p(R, R'; (3). However, the full density matrix of an in- 
teracting iV-particle quantum system is generally not known at low temperatures. 
Therefore one needs to resort to the discrete representation of the Feynman path 
integral formula for a low-temperature density matrix, which will be discussed in 
detail in Section ||. 

The finite-temperature nature of PIMC makes this a complementary approach 
to zero-temperature Monte Carlo methods, such as variational Monte Carlo, or 
Green's function-based methods such as diffusion Monte Carlo, which are reviewed 
in Chapter 1 of this volume, tl PIMC is currently the only numerical method capable 
of directly addressing finite-temperature superfluidity and the superfluid transition 
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in helium. In addition, unlike the zero-temperature methods, PIMC does not require 
the use of a trial function, requiring as input only the particle masses, numbers, and 
interaction potentials, in addition to the temperature and volume. Consequently, it 
is a numerically exact technique, and is independent of the trial function bias prob- 
lems that zero-temperature methods may suffer from. However, because PIMC pro- 
vides thermodynamic averages, state-specific information is generally not available. 
Although in principle the density matrix contains information on the full eigenspec- 
trum of the Hamiltonian H, extracting this requires the numerical inversion of a 
Laplace transform, au Such inversions are known to be notoriously difficult in the 
presence of Monte Carlo noise. Thus, to date, path integral Monte Carlo has pro- 
vided only very limited dynamical information of a direct nature. Nevertheless, it 
has provided critical microscopic input into dynamical models for physical systems 
in helium droplets, and in conjunction with zero-temperature, state-specific calcu- 
lations for these finite helium systems, PIMC has proven to be a powerful means of 
investigating the dynamic consequences of atomic scale structure of a superfluid. □ 

Feynman first applied the path integral approach to liquid helium in 1953, and 
provided a consistent clarification of the role of Bose permutation symmetry in the 
lambda transition of liquid helium, i In Feynman's original treatment, the multi- 
dimensional integrals of Eq. (0) were approximated analytically. In the late 1980's, 
Ceperley and Pollock subsequently devised a Monte Carlo scheme for the exact 
numerical evaluation of these multidimensional integrals, in which the combined 
configuration and permutation spaces were efficiently sampled using a multilevel 
Metropolis method. This allowed direct quantitative application of Feynman's 
path integral approach to the superfluid state of helium for the first time. 

Since then, the PIMC method has been applied to provide a quantitative descrip- 
tion of numerous bulk and finite bosonic systems. In addition to extensive studies 
of bulk helium, PIMC has now been employed in the study of ^He/^He mixtures, B 
of helium and molecular hydrogen droplets, Bu and of helium and hydrogen films on 
various surfaces. EJ^EJ In this work, we focus on the application of PIMC to quantum 
simulations of finite helium clusters, ^Hcat. Since we are primarily concerned with 
the bosonic isotope of helium, for the remainder of the chapter, we will denote "^He 
as simply He, unless explicitly stated. The field of helium cluster research has grown 
very rapidly over recent years due to new possibilities of inserting molecular probes 
and studying their properties. 113 An overview of the experimental work in this 
area is provided in Chapter 9 of this volume. Accompanying this rise in experimen- 
tal studies, there has been a correspondingly increased demand for complementary 
theoretical study of these finite quantum systems. 

The earliest PIMC simulation of pure Hejv droplets, made in 1989, demonstrated 
superfluid behavior for sizes N as small as 64. Ej This result, together with zero- 
temperature calculations of the size scaling for pure cluster excitation spectra made 
at that time, lla was taken to be strong theoretical evidence that these finite-sized 
clusters were indeed superfluid. This was later conflrmed experimentalhfthrough 
a series of elegant experiments with impurity-doped helium clusters. 1130 Doped 
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clusters present additional technical challenges for PIMC beyond the requirements 



C work with doped clusters ad- 
These studies showed that the 



posed by a finite cluster of pure HeAr. Early 
dressed the widely studied HcArSFg system, 
global superfluid fraction appeared not to be significantly modified by introduction 
of an impurity. However, it soon became apparent that interesting new local fea- 
tures due to Bose exchange symmetry were present in the immediate vicinity of 
an impurity. This led to the recognition that a local _non-superfluid density could 
be induced by the molecular interaction with helium. eI PIMC simulations of doped 
clusters have now been made with a variety of impurities, including OCS, HCN, ben- 
zene, H2, neutral and ionic alkali atoms, and some complexes of these molecules. As 
will be outlined here, these studies have revealed a broad range of properties of the 
dopant as well as insight into the dopant infiuencc on the superfluid properties of 
the droplets. Key features that have emerged from these studies of doped droplets 
are the ability to analyze superfluid behavior in nanoscale dimensions, to charac- 
terize quantum solvation in a superfluid, and to probe the atomic-scale behavior of 
a supcrfiuid near a molecular interface. 

In this review, we first provide an overview of the Feynman path integral formal- 
ism for quantum statistical thermodynamics in Section ^. Following this introduc- 
tion, we discuss the PIMC implementation of the general Feynman theory, focusing 
in particular on the multilevel Metropolis sampling method of Ceperley and Pol- 
lock. We then review applications of the PIMC method to the study of pure helium 
droplets, and of doped helium clusters containing atomic or molecular impurities in 
Section ^. The analysis of superfluidity in finite droplets, and concepts of nanoscale 
superfluids and local superfluidity are described in Section 3.4. PIMC applications 
to spectroscopic studies of doped helium clusters are summarized in Section ^ We 
conclude in Section ^ with a summary of open questions. 



2. Theory 

2.1. General formulation 

Here we deal with a cluster of N He atoms doped with a single impurity. In many of 
the studies made to date, the impurity is assumed to be fixed at the origin without 
either translational or rotational motion. Neglect of the impurity translational de- 
grees of freedom is not essential, but for impurity particles which are heavy relative 
to a helium atom, it is reasonable and often convenient to ignore the translational 
motion of the impurity. However, the neglect of the impurity rotational degrees of 
freedom should be treated with caution, especially for impurities with small prin- 
cipal moments of inertia. Incorporation of the rotational motion of the impurity is 
an area of current work. Hence for the present discussion we consider the following 
system Hamiltonian H: 

N N 
H^K+V^-Xj2\/^+Y^ 1/He-Hc(ry) + J2 ^Hc-imp(rO, (1) 
i—1 i<j i=X 
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where A = ?i^/2m4, with 7714 being the hehum mass. The potential energy V includes 
a sum of He-He pair potentials, Vhc-Hc, and He-impurity interactions VHc-imp, 
where the latter are most readily given in the molecular frame. If necessary, e.g. for 
light molecules, the impurity translational degrees of freedom can be incorporated 
by adding an additional term — A/Vj, corresponding to the impurity center-of-mass 
kinetic energy. 

In the path integral approach, one uses the identity e^^^'-^^^''^ — e^^^^ e^^'^^ 
to express the low-temperature density matrix by an integral over all possible paths, 
{i?, i?2, • ■ • , Rm-ItR'}, with the weight for each path given by the product of 
density matrices at a higher temperature T' = MT: 

p{R,R'-(3)^ j dRidR2 . . .dRM-i p{R,Ri;t)p{RuR2;t) . . . p{Rm-uR';t). (2) 

Here r = P/M = (kBT')~^ constitutes the imaginary time step defining the dis- 
crete representation of the path integral. For a sufficiently high temperature T' or, 
equivalently, for a small enough time step r, there exist several approximations to 
the density matrix that are sufficiently accurate for this factorization to be used in 
numerical work. The simplest of these high-temperature approximations is the 
primitive approximation, which is based upon the Trotter formula: Ell 

p(i?fe,i?fe+i;r) « J dR' {Rk\e-^^\R'){R'\e-^^\Rk+i). (3) 

The potential energy operator V in Eq. (|l|) is diagonal in the position representation, 
{R'\e--^\Rk+i) = e-^^f^'^+i^Jli?' - Rk+i), (4) 
while the kinetic term corresponds to the free particle density matrix: 

(i?,|e--^|i?') = p,{R,,R';r) = (47rAr) "^A'/^e" V4Ar^ (5) 

From Eqs. (^)-(^, the path integral representation for the density matrix in the 
primitive approximation may be expressed as: 

p(i?o,i?M;/3) = (6) 
dR^---dRM-iexp -EE -rY.V{Rk) 

; 1 — 1 J 1 



k=l i=l fc=l 



N N 



(4^Ar)-3A'W2 
with 

V{Rk) = E ^He-He(?'»j\fc) + E ^Ho-imp(r»,/c). (7) 

A single Rk is referred to as a time slice, and r^ fc, the position of the i^^ particle 
at the k^^ time slice, as a head. From here on, i will denote particle index and k 
will denote time index. Eq. may be viewed as a classical configuration inte- 
gral, with the exponent of its integrand corresponding to an energy function. The 
first term in the exponent, derived from the kinetic energy, corresponds to a spring 
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potential that connects beads representing the same atom at successive imaginary 
times, with coupHng constant (2Aidr^. This chain of beads connected by springs 
is often referred to as a polymer. □'c3 The hehum-hehum interaction is represented 
by an inter-polymer potential that has non-zero interactions only between beads 
located on different polymers and indexed by the same imaginary time value. This 
corresponds to Feynman's original idea of mapping path integrals of a quantum sys- 
tem onto interacting classical polymers, with a special form of polymer interaction 
potential. eI In the absence of kinetic contributions from the impurity, the helium- 
impurity interaction acts as an additional external field for the helium atoms, and 
hence for these polymers. 

The thermal density matrix p{R,R';(3) = {R\e-f^"\R') can be expanded in 
terms of the eigenvalues {£'«} and the corresponding eigenfunctions {0„} of the 
Hamiltonian H: 

p{R,R';(3) = ^0„(i?)C(i?')e-''''". (8) 

n 

This is appropriate for a system of distinguishable particles under Boltzmann statis- 
tics. For a Bose system such as He or para-H2, it should be symmetrized with respect 
to particle exchanges. This can be done by modifying the sum in Eq. (^) to a sum 
over exchange-symmetrized stationary states c/fa only: 

PB{R,R';P)=J2MRmR')e-^''''- (9) 

a 

A symmetrized eigenfunction (pa (i?) can be obtained by summing a stationary wave= 
function of distinguishable particles, 0„(7'i?), over all A^-particle permutations V: c3 

V 

Inserting Eq. jl^ ) into Eq. (^), we obtain the symmetrized density matrix for a 
Bose system: 

Pb{R, R'; /3) = ]^ E p(^^ '^)- (11) 
p 



2.2. Density matrix evaluation 

In order to make PIMC calculations more tractable, one wishes to use the smallest 
possible number of time slices M for a given temperature T. This means that it is 
essential to find accurate high-temperature density matrices at as small a value T' as 
possible, so that the imaginary time step r be kept as large as possible. It has been 
found that for the helium-helium interaction, the primitive approximation described 
in Eqs. is accurate enough at temperatures higher than ~1000 K. EJ This 

implies that a PIMC simulation at T ^ 0.3 — 0.4 K, where the spectroscooic mea- 
surements for the impurity-doped helium clusters have been performed, E3 would 
require several thousands of time slices. This computational expense can be avoided 
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by going beyond the primitive approximation to a more sophisticated approximation 
for the high-temperature density matrices. Based on the Feynman-Kac formula, E3 
the high-temperature density matrix p{R,R'\t) can be approximated by a prod- 
uct of the free particle propagator po{R,R';t) of Eq. and an interaction term 

g-(7(J?,ij';T). 

piR, R'- t) « po(i?, i?'; r)e-^(«'«'^-). (12) 

The interaction term e~'^(^'^ ''^^ is in turn factored into contributions deriving from 
the hclium-hclium and helium-impurity interactions, phc-Hc and PHc-imp, respec- 
tively. For spherical interactions, one can generate a pair-product form of the exact 
two-body density matrices using a matrix squaring approach discussed in detail in 
Ref. ^ We use this for the helium-helium interaction, which is spherically symmet- 
ric. Such helium-helium density matrices of the pairjjpmduct form have been shown 
to be accurate for /ks > 40 K, i.e., T' > 40 K. This same approach can be 
used for the helium-impurity interaction when this is also isotropic. However, for 
molecules the helium- impurity interaction VHc-imp(r) is in general not isotropic, 
and may involve complicated three-dimensional dependencies. This can be dealt 
with in several ways. One approach is to expand the helium-impurity interaction 
in spherical terms and then employ pair-product forms as above. We have found 
it convenient to work within the primitive approximation for the helium-impurity 
interaction, which allows considerable flexibility when changing impurities. The re- 
quired time step for the accurate primitive helium-impurity density matrices varies, 
depending on the impurity molecule involved {e.g., t^^ /ks > 80 K for He-SFg and 
He-OCS, /ks > 160 K for He-benzene). This must be recalibrated, e.g., by es- 
tablishing converged helium densities, for every new molecule that is studied. The 
same re-calibration requirement holds also for spherical expansions. 

2.3. Multilevel Metropolis algorithm 

For a diagonal operator O in the position representation, {R\0\R') = 0{R)S{R—R'), 
we need to consider only the diagonal density matrices for evaluation of its thermal 
average, Eq. (|l|). For the diagonal density matrix, both the sum over permutations 
in Eq. (^l]) and the multidimensional integration in Eq. can be evaluated by a 
sampling of discrete paths which end on a permutation of their starting positions, 
i.e., s = {Ro, Ri, R2, ■ • • , Rm-1, Rai} with Rm — VRq. This gives rise to an iso- 
morphic mapping onto ring polymers. In fact, all physical quantities discussed in 
this review can be estimated from a set of stochastically sampled ring polymers. 

In the sampling process, it will be natural to choose the probability density 
function as 

M-l 

= Z-^ W p{Rk,Rk+ur). (13) 

fc=0 

The Metropolis algorithm, a widely-used Monte Carlo sampling technique, provides 
a route to obtain the converged distribution tTs in the limit of many sampled config- 
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urations, as long as the detailed balance condition is satisfied for transitions between 
successive configurations: 

TT^Ps^s' = TTs'Ps'-*s- (14) 

Here Ps~^s' is the transition probability from a configuration s to s' . This is factor- 
ized into an a priori sampling distribution Tg^s' and an acceptance ratio As^s' ■ 

Ps^s' —Tg^giAg^gi. (15) 

In order to speed up convergence times in a path integral simulation, in partic- 
ular one involving permutation moves, it is very important to select an appropriate 
distribution function Tg^gi for a trial move s' from s. The most efficient way of 
doinfftbis is the multilevel Metropolis algorithm developed by Pollock and Ceper- 
ley. EjO Here one first chooses end points of each path by sampling a permutation 
v. Then the paths are bisected and the configurations at the midpoints sampled. 
This process of bisection and midpoint sampling is repeated multiple times, re- 
sulting in a multilevel scheme that samples whole sections of the paths in a single 
step. The acceptance ratio at each level of this multilevel Markov process is set so 
that the combined process of permutation and configuration moves may lead to the 



probability density function of Eq. (13). Detailed procedures are summarized as 
follows: c5 

(1) Initialize a configuration s. Typically one starts from a classical configuration, 
in which all beads representing each atom are located at the same site. So each 
polymer corresponds initially to a single point. 

(2) Choose a time slice k randomly between and M — I and construct a table 
for trial permutation transitions between time slices k and k + n, where n — 2^ 
and / is the level of this path updating process. For the simulation of a He 
system, I — 3 turns out to be a good choice for the permutation moves. Trial 
permutations may be restricted to cyclic permutations among 2, 3, or 4 particles. 
The probability for permutation transitions is proportional to 

-(i?fc — VRk+nY 



T-p = exp 



AXriT 



(16) 



Cj = J2Tv (17) 

V 

Thus the transition probability for permutational moves does not depend on the 
potential energy. Note that one can explore the entire iV-particle permutation 
space by repeatedly sampling cyclic permutations among a small number of 
particles. 

(3) Select a trial permutation V involving p atoms such that 

V'KV T"<V 

where x is a random number on (0, C/). This selects the permutation with 
probability T-p/C/. Then compute Ao = T-p/T/. After this, we will sample the 
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intermediate path coordinates connecting Rk with VRk+n- The coordinates of 
the {N ~p) atoms not on the cycle represented by V will not change from their 
old positions. This is level sampling. 
(4) Start a bisection algorithm by sampling a new midpoint R'l^j^n . For the sampling 
distribution function T{R'i._^n\Rk, VRk+n] jt): we use a multivariate Gaussian 
form centered at the mean position R — {Rk + Rk+n)/'^ (sec Eq. (5.16) of 
Rcf. ^). Then compute 

^ p{Rk,K+^;^r)p{R',^„, VRk+n; ^r) ^^^^ 

piRk, Rk+^; ^T)p{Rk+^, Rk+n] f t) 



Proceed to the next step with probability 

(20) 



AiT{Rk+R \Rk,Rk+n] f t) 



AoT(i?;,+ ^|i?fe,Pi?fc+„;fT)- 

If rejected, go back to step 3 and sample a new trial permutation. This is level 
1 sampling. 

(5) At the second level, sample R'k+n. and -Rj^^sn. by bisecting the two intervals 
and continue to the next level with the same procedures as used in step 4. This 
bisection process is repeated until we get to the final l-lh level. At the l-ih level, 
sample R'k+i-^ • ■ ■' ^'^'i R'k+n-i with the probability distribution function 

T{Rj\Rj_j^, Rj_^^;t). Proceed to the next step with probability 

A; ^ T{Rj\Rj_i,Rj+i:T) 



^'-1 T(i?^l^^-i,i?^+i;r) 



where 



A, . n (22) 

If rejected, go back to step 3. Fig. depicts the structure of a multi-level sam- 
pling with Z = 3, in a path integral containing M = 16 time slices. 

(6) Construct a new permutation table for all 2, 3, or 4 particle exchanges V' acting 
on VRk+n- Accept a new path Rk, R'k+ii • ■ • > R'k+n-ii'PRk+n with probability 
Ci/C-p, where C-p = J2v' Tvv- If rejected, go again back to step 3. 

(7) After replacing old coordinates and permutation table with new ones, repeat 
steps 3 to 6. 

(8) After attempting several hundreds or thousands of permutation moves between 
times slices k and k + n, followed by the bisection procedures to update the 
midpoints of all I levels, we select a new time slice k and repeat steps 2 to 6. 

One can check that the multilevel bisection algorithm described here satisfies 
the detailed balance condition in Eq. (^^. Note that Pg^s' is the total transition 
probability to go through all I levels. 
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2 



level 3 




k + n/4 



k + n/2 



levels 



k + n 



k + 3n/2 



levels 



\ 



level 3 



Fig. 1. Schematic of multilevel sampling. The figure shows a ring polymer of configuration beads 
for a single particle corresponding to M = 16 time slices, to be updated with a three-level (I = 3) 
sampling of 2"^ = 8 time slices simultaneously. The bold connections indicate the section of 8 time 
slices that is to be updated. 

2.4. Estimators for some physical quantities 

With the generahzed MetropoHs sampUng of the permutation symmetrized density 
matrix pb{R, R] (3), the thermal average of an observable O diagonal in the position 
representation can be estimated by taking an arithmetic average of 0{R) = {R\0\R) 
over the paths sampled. For instance, the helium density distribution around an 
impurity molecule can be estimated by 



Note that all time slices in ring polymers can be considered as equivalent. Unlike 
importance-sampled diffusion Monte Carlo methods, the PIMC calculation of struc- 
tural properties such as the density distribution does not involve any trial function 
bias. 

There are many ways to compute the energy in PIMC, discussed in detail in 
Ref. |2^. Most of the PIMC applications discussed here employ the direct estimator 
obtained by directly applying the Hamiltonian operator to the density matrix in the 
position space. For calculations neglecting the impurity translational and rotational 
degrees of freedom, the kinetic energy average is expressed in the path integral 
representation as 




A/-1 N 



(23) 



fc=0 i=l 




(24) 
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where = d/dRk is the SA^-dimensional gradient operator, and U'' = 

U{Rk-i, RkTT) is the interaction for link k, i.e., for the spring connecting beads 



k — 1 and k (see Eq. (12)). Computation of the potential energy is straightforward 



since this is diagonal in the position representation: 

- M E (^(^fc))- (25) 

fc=0 

As noted earlier, the inter-polymeric potential acts only between beads defined at 
the same imaginary time. 

One of the most interesting properties of bulk and finite He systems is their 
superfluid behavior. For bulk systems superfluid estimators are generally derived 
from linear response theory, i.e. by considering the helium response to boundary 
motion. Ej Pollock and Ceperley showed how to derive momentum density correla- 
tion functions that quantify the superfluid response of bulk systems with periodic 
boundary conditions. Sindzingrc et al. subsequently developed a global linear re- 
sponse estimator for finite helium clusters with free boundaries. 113 This estimator 
is based on the response to a rotation of continuous angular frequency, i.e. to a 
classical rotation such as might be appropriate to a macroscopic droplet. Consider 
the Hamiltonian in a coordinate frame rotating about an axis with frequency to, 

Hrot =Ho-t-uJ, (26) 

where L is the total angular momentum operator. For a classical fluid, in the limit of 
an infinitesimally small rotation the entire fluid should rotate rigidly with classical 
moment of inertia Id. But in a Bose superfluid, only the normal component responds 
to the rotation, resulting in an effective moment of inertia 



d{t) 



duj 



(27) 



Note that this is to be evaluated in the limit ^ 0, appropriate to rotation of a 
macroscopic system. In a homogeneous system the normal fraction can be defined 
as 

^ = (28) 

P J^cl 

and the complementary superfiuid fraction is then given by 

Pl = l-P^ = hL^_ (29) 
p p Id 

For uj = ujz, this linear response estimator can be expressed in the path integral 
representation as lil 

^ Am' {AD 

p ph'ici ' 
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with the vector quantity A defined as 

A = ^^(r,,fc X r,.fe+i). (31) 

i,k 

Here the summation runs over particle i and imaginary-time slice k. The vector 
quantity A is the directed total area of closed imaginary-time polymers spanned by 
all N particles, e.g.^ is the projection of A on the i-direction. The average size 
of a single polymer is given by its thermal de Broglie wavelength A = {[3h? /ra)'^/'^ . 
This becomes negligible in the high-temperature limit, and thus the corresponding 
superfluid fraction Ps/ p goes to zero (although the projected area can remain finite 
at the microscopic level). At low temperatures, when the de Broglie wavelength 
becomes comparable to the inter-polymer spacing, particle exchanges cause poly- 
mers to cross-link and form larger ring polymers. The projected area A increases 
correspondingly and the helium system attains an appreciable superfluid fraction. 

We discuss the application of this finite-system superfluid estimator in detail in 
Sec. ^, together with analysis ofhaw-f his can be decomposed into local contributions 
for an inhomogeneous system. BllaO 



3. Superfluidity and quantum solvation of atoms and molecules in 
bosonic helium clusters 

Spectroscopic studies of impurity-doped clusters have^Uowed experimental inves- 
tigation of a variety of excitations in helium clusters. Ej The relevant temperature 
range currently accessible is T ~ 0.15 — 0.5 K. Thus, the incorporation of Bose 
symmetry is essential in simulation of these systems. In this section, we focus on 
the application of finite-temperature PIMC to bosonic helium clusters. We begin 



by briefly reviewing studies of pure clusters in Sec. 3.1, and then focus on the more 



recent work for the clusters doped with various impurities in Sees. 3.2-3.4. Sees. 3.2 



3.3| summarize the structural and energetic aspects, while Sec. deals with the 
microscopic analysis of superfluid properties of the doped clusters. 



3.1. Pure clusters 

Most of the previous theoretical studies involving pure clusters are based on zero- 
temperature methods, and have focused on the cluster elementary excitation spec- 
trum, which qualitatively retains the phonon-roton features characteristic of bulk 
He H. Current w ork in this area aims to understand the physical nature of the 
roton excitations. Ea Zero- and finite-temperature calculations for pure helium clus- 
ters have been reviewed previously, and so we shall provide only a brief outline 
of the finite-temperature results here. The first studies were made by Cleveland 
et al, usinff the path integral molecular dynamics approach in which exchange is 
neglected. t2l This allowed structural analysis, which was used to study the changes 
in droplet density and diffuseness as a function of size. Permutation exchange sym- 
metry was incorporated by Sindzingre et al. in a study of the temperature and size 
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dependence of the global superfluid fraction in finite Hcat clusters. These calcu- 



lations employed the area estimator discussed in Sec. 2.4 and showed that a broad 
transition to a predominantly superfluid state occurs at a temperature depressed 
from the bulk superfluid transition temperature, in accordance with expectations 
from scaling of phase transitions for flnite systems. The extent of depression in- 
creased as the cluster size decreased. For N — 64, the onset of the transition occurs 
just below r = 2 K, and the transition appears complete at T = 0.5 K, with about 
90% or more of the cluster being superfluid at that temperature. A qualitative ex- 
amination of the relative contribution of long exchange path lengths to the density 
revealed that the long exchange path contribution was largest in the interior of the 
droplets. 



3.2. Atomic impurities 

Neutral and ionic atomic impurities constitute the simplest dopants. For ground 
electronic states, the helium-impurity interatomic potential can be calculated with 
fairly high accuracy using standard quantum chemistry methods, and the helium- 
helium interatomic potential is well-known. Thus, within the two-body approxima- 
tion, it is possible to construct accurate potential energy surfaces for the ground 
electronic state. The interactions of excited electronic states with helium are, by 
comparison, less well-characterized and only a few calculations of electronically 
excited potential energy surfaces have been even attempted. To date, PIMC calcu- 
lations have been made for the neutral alkali metal impurities Li, Na, K, and for 
the ionic impurity Na^ . E3 In general, the solvation characteristics p£ each impurity 
are controlled by a balance between diff'erent energetic factors. Ej'Ea These include 
the helium-impurity interaction strength, the helium-helium interaction strength, 
the impurity kinetic energy (and thus impurity mass), and the free energy change 
due to the loss of exchange energy for helium atoms adjacent to the impurity. 

The He-Li, _He-Na, and He-K ground state potentials typically_have well depths 
of ^ 1 — 2 K, E3 smaller than the He-He well depth of ~ 11 K. 123 By considering 
these potential energy factors alone, one would qualitatively expect that the atomic 
impurities would reside on the droplet surface in order to minimize the total energy. 
The PIMC studies, made at T = 0.5 K, indicate that the neutral alkali impurity 
species are indeed surface-attached for cluster sizes of iV < 300. EJ In all cases, the 
perturbation on the cluster structure due to the presence of the impurity is weak. 
The neutral impurity atom induces small but distinct modulations in the helium 
density, starting at the surface and decaying into the interior of the droplet. While 
no zero-temperature microscopic calculations of these systems have been made to 
date, it is expected that this behavior would persist to lower temperatures and is 
therefore also applicable to the experimental studies of these alkali-doped systems 
that are made at T = 0.38 K. @ 

The ionic impurity Na^ interacts much more strongly with helium and conse- 
quently gives a markedly stronger structural perturbation of the local helium den- 
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Fig. 2. Helium density profiles relative to the impurity center, at T = 0.625 K. The left panel (a) 
shows the radial density profiles for HeiooNa^ (dashed lines) and He64SF6 (solid lines), where an 
isotropic He-SFg interaction was used. The right panel (b) shows anisotropic ifeelium densities for 
SFe (solid lines) from a calculation using an anisotropic He-SFg interaction, viewed along the 
molecular C3 (higher values) and C4 (lower values) axes— This is compared to anisotropic radial 
density profiles for Heag-benzene (dotted-dashed lines), viewed along the molecular Ce axis 
(higher peak) and along the C— H bond (lower peak). The HeiooNa+ profile is reproduced from 
Ref. ||. 



sity. The He-Na+ well depth is 407 K, about 40 times larger than the He-He well 
depth. Here the finite-temperature PIMC studies indicate that the Na"*" ion resides 
in the center of the cluster, and the strength and range of tke He-Na"*" interaction 
induces a tightly packed helium "snowball" around the ion. E3 In Fig. ^a, the radial 
helium density profile for the HciooNa"'" system at T = 0.625 K is shown. This 
is compared with the iV = 64 radial density profile around the molecular impurity 
SFg, at the same temperature using an isotropic He-SFe interaction. There is a very 
strongly modulated layer structure around the Na^ ion, with a high first coordina- 
tion shell peak followed by a second peak of lower density. Similar structural features 
have b e e n seen in variational shadow function calculations for Na"*" and K"*" in bulk 
He, although quantitative differences exist in comparison with those results. 
In the variational calculations the local angular ordering within the coordination 
shells was also examined, leading to more conclusive evidence of solid-like structure 
in the first two shells. These studies indicate that there definitely exists a more 
strongly layered shell structure in the helium density around an impurity ion than 
around neutral atomic species, with more solid-like character. This feature can be 
further explored in the imaginary-time path integral representation by examining 
the permutation exchanges of helium atoms at specific locations. For the HciooNa+ 
cluster at r = 0.8 and 1.25 K, the atoms in the first coordination shell rarely partic- 
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ipated in permutations with other particles, and thus are well-localized in the PIMC 
sense. In the second solvation shell, some atoms are involved in long exchanges at 
the lower temperature, while in the outer third shell most atoms are involved in 
long exchanges. From this, it was inferred that the third shell is superfluid, while 
the second shell has an intermediate, temperature-dependent character. E3 Such an 
analysis has been made in more quantitative detail for the molecular impurities, 
which we discuss next. 

3.3. Molecular impurities 

Molecular impurities introduce an additional level of complexity because molecules 
have internal structure and usually possess an anisotropic interaction with helium. 
Especially for the larger molecules, there is a severe lack of accurate two-body 
molecule-helium interaction potentials. Nevertheless, the study of molecular impu- 
rities in helium clusters is currently of great interest, with an increasing number of 
experiments being performed on a variety of molecules. Even with simple models 
for the molecule-helium interaction, analysis of these experiments in terms of the 
perturbation of the helium environment on the molecular internal degrees of free- 
dom has provided much insight into the quantum fluid nature of these clusters. It 
is important to recall that to date, all PIMC work involving molecular impurities in 
helium have not explicitly incorporated the impurity rotational kinetic energy. This 
is not an essential restriction, and has been so far made for convenience rather than 
for any fundamental limitations. Since zero-temperature DMC calculations have 
recently shown that the helium densities ammicLsmall molecules may be sensitive 
to the rotational motion of the molecule, bM^eS it would be desirable to incor- 
porate the rotational degrees of freedom in future path integral studies. Only for 
the heaviest rotors can molecular rotation be justifiably omitted. Several studies 
have also justified neglecting the translational motion of the molecular impurity, in 
which case the helium atoms may be regarded as moving in the external potential 
field of themolecule in the molecular body-fixed frame, given by the Hamiltonian of 
Eq. (^. eH22I The validity of this assumption can be assessed with the comparative 
study of molecular dfilpcalization as a function of molecular mass and binding en- 
ergy shown in Fig. ^. CJ There, the probability of finding a molecule at some distance 
r from the cluster center-of-mass is shown for a series of molecules. As expected, 
the heaviest dopants such as SFg and OCS arc well-localized near the center of the 
cluster, and thus it is a reasonable approximation to neglect their translational mo- 
tion. On the other hand, H2 is much more delocalized throughout the interior of the 
cluster, due to its comparatively smaller mass and weaker helium-impurity interac- 
tion, and therefore requires that its translational motion be properly incorporated. 

To date, the most extensively studied molecular impurity is the octahedral SFg 
molecule. Early PIMC work on-SFg in helium clusters employed isotropic molecule- 
helium interaction potentials, 113 and was later extended to include anisotropic in- 
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Fig. 3. Radial distribution for several impurity molecules (H2, HCN, SFg, and OCS) relative to 
the cluster center-of-mass, shown as P{r) = 47rr^p(r) such that J P(r)dr = 1. All calculations 
were made from PIMC at T = 0.312 K, and include the impurity translational kinetic energy. 
Isotropic interaction potentials were used, and Bose permutation symmetry was not included. The 
H2 and HCN distributions were obtained from a calculation with N = 128 He, while the SFq and 
OCS distributions correspond to a cluster of TV = 100 He. Data courtesy of D. T. Moore. 

teractions. The heUum structure around SFg from an isotropic calculation is 
shown in Figs. ||a and_J|. The anisotropic He-SFg potential surface has a global 
minimum of —84 K, c3 considerably deeper than that of the He-He interaction. 
Thus, SFg is expected to reside aiJihe center of the cluster. This has been verified 
by both zero-temperature DMC e3 and by finite-temperature PIMC calculations. 
For the He64SF6 cluster in the temperature range of T = 0.3 — 0.75 K, there is an 
anisotropic layering of the helium density around the SFg. Integration of the helium 
density over the first solvation shell yields about 23 atoms, independent of whether 
isotropic or anisotropic interactions are employed. □ The strength and range of the 
molecule-helium interaction pins the helium density in the first solvation shell to 
a total density comparable to that of the more strongly bound HeTv-Na"*" system. 
Detailed analysis of the helium density distribution around the molecule shows that 
while the angular average of the density in the first solvation shell is independent of 
temperature below T = 1.25 K, there is a small increase in the extent of anisotropy 
as the temperature is lowered. This is illustrated in Fig. ^, with a comparison of 
the densities along different molecular symmetry axes for an = 64 cluster at 
temperatures T = 0.625 K and T — 0.312 K. As the temperature is increased above 
T = 1.25 K, this trend to an increasingly isotropic distribution is further modified 
by the onset of evaporation of helium atoms. Evaporation begins with atoms in the 
second solvation shell, is clearly evident at T — 2.5 K, and is essentially complete 
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at T 5.0 K (Fig. Isl) 
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Fig. 4. Comparison of the helium density distribution around the octahedral SFe molecule in 
a Af = 64 cluster at T = 0.625 K and T = 0.312 K. Solid lines show the lower temperature 
densities and dashed lines the higher temperature densities. Panel (a) shows the angular- averaged 
density poi for which the profiles at two different temperatures are identical. Panels (b), (c), and 
(d) show the densities along the three symmetry axes of the molecule C2 , C3 , and C4 , respectively. 
The higher temperature profiles show consistently smaller peak values in the first solvation shell, 
indicating a decrease in the anisotropy of the distribution as temperature increases. 



Since the first experiments for doped clusters that employed SFg as a probe 
species, @ a hipad array of molecular impurities have been studied by spectro- 
scopic means. l3 The infrared spectral regime has provided a particularly rich field 
of study. Vibrational spectra in the infrared at T ~ 0.4 K show rotational fine 
structure in *He droplets, but not in '^He droplets, providing evidence that quan- 
tum statistics play an important role in the spectral properties of the dopant. There 
is now an increasing collection of experimental data available for the rotational dy- 
namics of molecules possessing varying symmetries and a range of values for the gas 
phase rotational constant. To date, PIMC has been used to make theoretical stud- 
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Fig. 5. Helium evaporation for the SFsHesg cluster. The lower panel plots a snapshot of imaginary 
time paths at T = 1.25 K. At this temperature the helium atoms arc bound to the cluster. In the 
middle panel, at T = 2.5 K, the cluster begins to dissociate, loosing helium atoms. In the upper 
panel, at T = 5.0 K, the cluster has completely evaporated. 
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ies of thejjnear rotors OCS i@ and HCM-i the planar aromatic molecule lienzene 
(CeHe), S the linear (HCN)3 complex, and the 0CS-(H2)m complex. From 
these studies the notion of two different dynamical regimes has emerged, namely 
that of heavy molecules such as SFg that are characterized by gas phase rotational 
constants Bq < 0.5 cm~^, and a complementary regime of lighter molecules pos- 
sessing larger gas phase values of Bq. a This division into two dynamical regimes 
based on rotational constants emerges from analysis of the helium solvation density 
and energetics derived from path integral calculations. 

The OCS impurity lies in the regime of relatively heavy molecules, with Bq = 
0.20 cm~^. The He-OCS potential has a global minimum of ~ 64 K, c3 which is only 
about two-thirds that of the He-SFg potential, o It is important to consider the 
anisotropy of the intermolecular potential in addition to its strength when assessing 
the quantum solvation structure. In this respect the linear OCS molecule has lower 
symmetry than the octahedral SFg. The minimum angular barrier for rotation of 
the OCS about an axis perpendicular to the molecular axis {i.e., the angular adi- 
abatic barrier for rotation) is 41.9 K. @ This barrier is markedly higher than the 
corresponding value 20.7 K for SFg, and consequently gives rise to stronger angular 
modulations in the solvating density, i As shown in Fig. ^, PIMC calculations for 
the IIe640CS cluster at T = 0.312 K reveal a strongly structured helium density, 
forming approximately elliptical solvation shells around the OCS impurity. The first 
shell integrates to ~ 17 atoms, i Because of the axial symmetry of the Hc-OCS po- 
tential, the density at the global minimum forms a ring around the OCS molecular 
axis, consisting of about 6 helium atoms. 

The benzene molecule (CgHg) also lies in the heavy regime. The benzene tt- 
electron character leads to a highly anisotropic interaction with helium, with two 
deep, equivalent global potential minimum locat ed nn the six-fold axis of sym- 
metry above and below the plane of the molecule. EZlc3 A PIMC study of benzene- 
doped clusters has shown a highly anisotropicJielium structure around the impurity 
molecule that reflects this six-fold symmetry. 123 The sharpest density peak is located 
along the Cg-axis, at the two equivalent locations of the global potential minima. 
These two global density maxima are higher than the local density maxima viewed 
along the in-plane directions by more than a factor of four, reflecting the marked 
anisotropy of the He-benzene interaction potential. The extreme density anisotropy 
is summarized in Fig. |^b where the dotted-dashed lines show density profiles along 
the Cg-axis and along one of the in-plane directions. Integration over any one of 
the two equivalent global density maxima gives exactly one helium atom. We see an 
interesting effect of near complete localization of these two helium atoms located 
at the two global potential minima on either side of the molecular plane. As noted 
in Ref. this phenomenon can be viewed as a precursor form of helium adsorp- 
tion onto a molecular nanosubstrate. Extending these studies to larger polyaromatic 
molecules will allow contact to be made with PIMC studies of helium adsorption 
on graphite. Il3 

In contrast to this highly structured quantum solvation observed around the 
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Fig. 6. Total helium dtmsity around HCN (top panel) and OCS (bottom panel) for a = 64 
cluster at T = 0.312 K. □ The origin is set at the impurity center-of-mass. The OCS is oriented 
with the oxygen end directed towards the -z direction, and the HCN is oriented with the nitrogen 
end directed towards the -z direction. 

heavier molecules such as OCS and benzene, the linear HCN molecule falls into 
the light molecule regime, with a significantly larger gas phase rotational constant. 
For HCN, Bq = 1.48 cm"i. The He-HCN potential E3 is both weaker (its global 
minimum is —42 K) and less anisotropic than the He-OCS potential. While there 
is clearly still an ellipsoidal layering of the helium density around the HCN, within 
each solvation shell there is now a noticeable lack of angular structure, in contrast 
to the situation with OCS (Fig. ^). For such a light rotor, neglect of the molecular 
rotational kinetic energy now becomes a more serious concern. From DMC stud- 
ies assessing the effect of molecular rotation, Bl30 the expectation here is that 
the helium density will become more diffuse when molecular rotation is explicitly 
incorporated into PIMC. 

Self-assembled linear chains of polymeric (HCN)jv/ have been detected experi- 
mentally in helium droplets. El The helium structure around such linear chains has 
recently been addressed with a study of the properties of helium droplets with up 
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to = 500 atoms that contain (HCN)2 dimers and (HCN)3 trimers. Like the 
monomeric molecules discussed above, the HCN polymers are found to be located at 
the center of the droplet and to induce a layering of the helium density. Draeger et 
al. have analyzed the structure of the first solvation shell around the linear polymer 
in terms of a two-dimensional film, estimated the effective confinement potential 
for displacement away from ihe droplet center, and made calculations for vortex 
formation in these droplets. It has been suggested earlier that the presence of 
a linear inmurity species might stabilize the formation of a vortex line in helium 
droplets. c3 The expectation here is that a vortex line could be pinned along the 
molecular axis of a linear molecule such as HCN, or more likely, along the axis of 
a linear polymeric chain such as (HCN) a/. While the physics of vortices constitute 
an essential feature of bulk He II, Ell and ways of producing and detecting vor- 
tices during helium droplet formation have been the subject of much discussion (see 



Ref. |53| and therein), no experimental evidence has been found so far for existence 
of vortices in finite helium droplets. Theoretically, vortices have been shown to be 
unstable in pure droplets, El and the situation with regard to doped droplets is still 
controversial. The energy for formation of a vortex, AEy, is defined as 

AEv =Ev-Eo, (1) 

where Eq is the ground state cluster energy and Ey is the energy of the cluster with 
a vortex line present. Within the fixed-phase approximation, the PIMC estimate for 
this vortex formation energy is ~ 30 K for a He5oo(HCN)M cluster aX T — 0.38 K, 
where M = — 3. In this case, the vortex formation energy is found not to be 
significantly affected by the presence of a linear impurity. In comparison, density 
functional calculations made for a range of impurities and cluster sizes give values of 
AEv that are larger than the fixed-phase PIMC estimates by^ factor of 3, and that 
are reduced by ~ 5 — 10 K in the presence of an impurity. An exact estimator 
for the energy of a cluster in an angular momentum state m relative to m = 
has been derived using angular momentum projection methods, cll Application of 
this estimator at T = 2.0 K indicates that the presence of an impurity actually 
results in a slight increase in the vortex formation energy. More work is required 
in this direction, in particular the systematic examination of the cluster size and 
temperature dependence of AEy obtained from the angular momentum projection 
estimator. 

Many other complexes have now been synthesized in helium droplets. 113 Indeed, 
these droplets are proving to be a remarkably versatile quantum matrix environment 
for synthesis of unusual or metastable aggregates. Of particular interest from a 
fundamental point of view are the complexes of OCS with molecular hydrogen, H2. 
Recent spectroscopic measurements on 0CS(H2)j\/ complexes inside Hcat clusters 
have shown an interesting feature that has been interpreted as evidence of nanoscale 
hydrogen superfluidity. Ell Initial PIMC studies of these systems have been carried 
out cJ using accurate pair potentials of OCS with He and with H2. oEj Since the 
H2-OCS potential surface has a similar angular modulation as that for Hc-OCS, 
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but a deeper minimum, the OCS molecule is expected to bind prefepaitially to H2 
over He. Calculations for the OCS(H2)5 complex in the Heag cluster eZI showed that 
approximately six helium atoms, which would normally occupy the region of the 
global potential minimum in the absence of i72, are completely displaced by five 
H2 molecules. These H2 molecules form a complete ring encircling the linear OCS 
molecule at the region of lowest potential energy. The helium density is pushed 
either to either the secondary peaks in the first shell, or outwards from first to 
second shell region. 

3.4. Exchange permutation analysis and impurity-induced 
non- superfluidity 

In addition to providing structural and energetic information, PIMC is currently the 
only numerical method capable of providing information on finite-temperature su- 
perfluidity in He systems. At high temperatures an A^-body system may be described 
by Boltzmann statistics, i.e. in the path integral representation, only the identity 
permutation is important. At low temperatures however, permutations must be 
included in the path integral representation for the thermal density matrix. In par- 
ticular, for liquid helium near the lambda transition, Feynman qualitatively showed 
that the presence of long exchange cycles gives rise to the sharp increase in the heat 
capacity, but due to the analytical approximations made in his analysis he was not 
able to correctly identify the order of the transition, i Further refinements in this 
and numerical PIMC simulations have quantitatively confirmed both the transition 
temperature and its order. c3 

The area estimator of Eq. ( |30| ) gives a scalar value for the global superfluid frac- 
tion Ps/ p. This provides a complete description for homogeneous helium systems. 
However, a finite cluster of nanoscale dimensions necessarily contains inhomogene- 
ity deriving from the surface, and atomic and molecular dopants provide additional 
sources of inhomogeneities. In this situation Eq. ( |30| ) may be interpreted as provid- 
ing an estimate of the global superfluid fraction averaged over all sources of inhomo- 
geneity. It is notable that the impurity molecule does not significantly perturb this 
global superfluid fraction. For neutral Na^^ioped clusters, the area estimator yields a 
global superfluid fraction of about 95%, cil consistent with the very weak perturba- 
tion of the density noted earlier. For more strongly bound systems such as HejvSFg 
and HcAfNa''", it is found that ps is similarly large, approaching unity for A^ > 100 
at r = 0.625 K. HH Thus to see a molecular effect on superfluidity, one needs to 
examine the local solvation structure on microscopic length scales. Here, the density 
p is no longer uniform, particularly in the neighborhood of an inhomogeneity. Thus, 
the superfluid fraction Ps/ p is expected to be dependent on position. Some indirect 
indications of this have also been found in analyses of helium films. 0^0 

A simple way to qualitatively estimate the local dependence of superfluid char- 
acter is to examine the probability np(r) of a particle at a position r to participate 
in an imaginary-time exchange cycle of length p. As discussed previously, Bose su- 
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perfluidity is associated with the existence of exchange cycles of long p. In a pure 
cluster, the single source of inhomogeneity is the cluster surface. For a pure clus- 
ter, np>6(r) goes to zero as the radial distance r approaches the surface. In the 



presence of an impurity, the examples discussed in Sees. p.2| and B.3 show that 
an embedded molecule can significantly modify the total density distribution p{r). 
Consequently one also expects changes in the local superfiuid character. Kwon and 
Whaley have systematically examined np(r) for helium clusters doped with a single 
SFg impurity. 113 They define a local superfiuid density by 

N 

Ps{r) ^ np(r)p(r) (2) 

p>p' 

where p(r) is the total density at r, and p' is a cutoff value for the permutation 
cycle length. This does not account for the tensor nature of the superfiuid response, 
providing a three-dimensional anisotropic representation of a scalar, that may be 
viewed as an average over the set of tensorial response functions. For a molecule 
with high symmetry such as a spherical top, this will not be a serious limitation. 
For linear molecules it will introduce some uncertainty. For clusters of > 50, 
most of the polymers sampled involve either one or two atoms {p = 1,2) or many 
atoms (large p). Thus in this size regime a clear cutoff exists. For these sizes, Kwon 
and Whaley used a value of p' = 6. For small clusters {N < 50), a clear distinction 
between short and long exchange cycles cannot be made, which implies that in the 
small cluster regime a two-fiuid interpretation of the density cannot be applied. 

For the octahedral SFg molecule, the local superfiuid estimator of Eq. (||) yields 
an anisotropic superfiuid solvation structure around the impurity molecule, whose 
density modulations are similar to those of the total density p{r). The local non- 
superfiuid density Pn(j) — p(r) — ps{r) does depends weakly on temperature, which 
implies that /?„ (r) consists of thermal contribution and a molecule- induced compo- 
nent. Fig. shows a three-dimensional representation of the local non-superfiuid 
distribution around the octahedral SFg molecule. The red areas of highest non- 
superfiuid density are located at the octahedral sites of strongest binding to the 
molecule, reflecting the origin of this as a molecule-induced non-superfiuid. This is 
in contrast to the thermal normal densily of bulk He II in the Landau two-fluid de- 
scription of a homogeneous superfluid. Ej The molecule-induced density component 
depends on the strength and range of the helium-impurity interaction potential, 
and is expected to persist at T = 0. Detailed analysis shows that it is non-zero only 
in the first solvation layer around the molecule. Ilj 

An analysis using the local estimator of Eq. ^ has been applied to a number 
of different molecular impurities in helium clusters, including the linear molecules 
OCS and HCN, eIeZI and benzene. EZl These systems exhibit a similar layering in 
both local superfiuid density ps (r) and local non-superfiuid density pn (r) around the 
molecule. The non-superfiuid density shows slightly stronger modulations, result- 
ing in a weakly anisotropic local superfiuid fraction in addition to the component 
densities themselves. eI In the more strongly bound He-OCS case, the maximum 
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Fig. 7. Local non-superfluid density pn{r) around SFe in a W" = 64 cluster at T = 0.312 K, as 
measured by the exchange path decomposition of the density. B The color scale goes from red for 
highest values of pn(r), to blue for the lowest values of Pn(r). The size of the ball corresponds to 
a distance from the SFe molecule of r = 9.0 A. The two cuts display the density in two equivalent 
planes containing C3 and C2 axes. The strong binding to the octahedral sites located along the 
C3 axes is evident, with 4 of the 8 octahedral sites visible here. 

of the non-superfluid component is roughly ~ 50% of the total density, while for 
the weakly bound He-HCN, the non-superfluid, or short-exchange path, component 
is only ~ 20%. We note that the molecule-induced non-superfluid density is also 
present around an impurity possessing an isotropic interaction with helium, i.e., it is 
not essential to have an anisotropic interaction. In fact the existence of a molecule- 
induced non-superfluid density was first seen in calculations of the SFg molecule 
with isotropic interactions potentials, summarized in Fig. ^. 

Nakayama and Yamashita have pursued a similar analysis of the local superfluid 
densiiy for the HeArNa"*" cluster, which exhibits a triple-layer structure for N ~ 
100. e3 While they did not explicitly compute the local quantities Ps{r) or p,i(r) in 
their PIMC study, they observed that the helium atoms in the first solvation shell 
(r < 4 A) rarely participate in long exchanges. This observation, combined with 
the pair distribution fmictions computed with respect to atoms in the first shell, 
led them to conclude that the first shell is solid-like. 
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Fig. 8. Total, local non-superfluid, and local supcrfluid densities around SFe in a Af = 64 cluster 
at T = 0.625 K, calculated with only the isotropic component of the SFe-He interaction potential. 
The origin is set at the impurity center-of-mass. The local superfluid density is calculated with 
the exchange path length criterion of Ref. ^| 

As discussed previously, an even more anisotropic impurity-helium interaction 
potential is provided by the benzene molecule. For the Heag-benzene cluster at 
T = 0.625 K, the two atoms corresponding to the two total density maxima localized 
at the two global potential minima, undergo less than 2% permutation exchanges 
with the surrounding helium. This implies that they are effectively rernQved from 
the superfluid, i.e., constitute a true "dead" adsorbed pair of atoms. Ej^eII This 
near-complete removal of individual helium atoms in the solvation shell from par- 
ticipation in permutation exchanges of nearby helium atoms has not been seen for 
other molecules to date. It provides an extreme case of the local non-superfluid den- 
sity, Pn(r), in which there is no longer any partial exchange of helium atoms between 
the non-superfluid and superfluid densities. These features of the helium solvation 
around a benzene molecule are expected to appear also in clusters containing larger 
polyaromatic molecules such as tetracene and naphthalene. A systematic analysis 
of the effect in planar aromatic molecules of increasing size, making the transition 
from a molecular to a micron-scale substrate, would be very useful. 

Recently another local estimator of superfluidity has been proposed that de- 
composes the projected area into contributions from each local density bin. E3 This 
decomposition allows the anisotropy of the response tensor to be evaluated ex- 
plicitly. Application of this local estimator to the linear HCN trimer embedded in 
helium droplets has confirmed that the superfluid density is reduced in the first sol- 
vation layer, consistent with the presence of a local non-superfluid density induced 
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by the molecule-helium interaction, as first established by Kwon and Whaley. El 
Furthermore, this new estimator shows that there is an asymmetry between the 
helium response to rotation about the molecular axis, versus rotation about an axis 
perpendicular to the molecular axis. Draeger et al. find that the superfluid response 
is reduced more foiixotation about the perpendicular axes than for rotation about 
the molecular axis. EZl In both cases it is less than unity, implying that there is a non- 
superfluid component both when rotation is accompanied by variation in potential 
energy, and when there is no variation in potential energy. This finding supports 
the existence of a local non-superfluid induced by an isotropic helium-impurity in- 
teraction, using the exchange path analysis of Kwon and Whaley (Fig. ||). Thus 
the local non-superfluid is not dependent on the presence of anisotropy, but derives 
primarily from the stronger attraction of helium to the molecule than to itself. 

These studies of various molecules embedded in He at clusters employing differ- 
ent estimators of local superfluidity all point to the existence of a molecule-induced 
non-superfluid density in the first solvation shell around a molecule. While the de- 
tails of this non-supcrfluid density may be somewhat dependent on how it is defined, 
it is evident from the studies of OCS, benzene, and HCN polymers made to date, 
that this local non-superfluid component is a general phenomenon to be expected 
for all heavy molecules. It therefore appears to be one of the defining features of 
quantum solvation in a superfluid. The extent of exchange between non-superfluid 
and superfluid densities exhibits a dependence on the strength of the helium inter- 
action with the molecule. Benzene provides an interesting extreme case of negligi- 
ble exchange between non-superfluid and superfluid density components, while less 
anisotropic molecules such as SFg still possess considerable exchange between local 
non-superfluid and local superfluid. Thus, both the interaction strength with the 
molecular impurity and the symmetry of this interaction are important. The ben- 
zene example indicates that there are useful analogies with the well-known "dead" 
or "inert" layer of helium adsorbed into bulk solid surfaces, which will be valuable 
to pursue in future studies. 



4. PIMC and the connection to cluster spectroscopy 
4.1. Electronic spectra in J?ejv 

Calculations of electronic spectra typically require accurate potential energy sur- 
faces for both ground and excited electronic states. This is particularly challenging 
for excitations in condensed phases. To date, theoretical work in this area has been 
limited to relatively to simple systems, where the helium-impurity ground and ex- 
cited state pair potentials can be computed to good accuracy using standard ab 
initio electronic structure methods. Thermally-averaged electronic absorption spec- 
tra for the ^£ <— '^S transition have been computed for neutral alkali impurities at 
T — 0.5 K, c3 using a modification of the semi-classical Frank-Condon expression 
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for the electronic lineshape, 

I{u;)(x\M\^ J dRp{R,R;(3)S[Ve{R)-Eg{R)-hu;], (1) 

where M is the electronic transition dipolc moment, and Ve is the potential in the 
electronic excited state. The quantity Eg is a local ground state energy, which is 
assumed to take the form 

N N 

Eg{R) = Ti„p(i?) + ^ FHc-imp(r«) + J2 ^Ho-Hc(r.y ), (2) 

2—1 i<j 

and explicitly incorporates the kinetic energy of the impurity atom Tjmp. The terms 
l^e-imp and Vhc-Ho correspond to the helium-impurity and helium-helium ground 
state pair potentials, respectively. The electronic excited state potential Ve is ob- 
tained from the diatomics-in-molecules (DIM) model, 

N 

Ve{R) = F4-imp(^) + E ^He-He(r.,), (3) 

i<j 

where the first term V^^_^^^^p is the adiabatic energies of the alkali atom in the 
manifold interacting with the N helium atoms, and the remaining helium-helium 
pair potentials Vhc-Ho are taken to be identical to the ground state. Thus, the 
thermal absorption profile can be computed by sampling this energy difference 
of Eq. (|^) from a PIMC simulation. As discussed in Sec. 3.2, the neutral alkali 



impurities reside on the droplet surface, and the resulting perturbation on droplet 
properties is weak. The electronic lineshape is therefore most sensitive to the details 
of the surface structure near the alkali atom. The PIMC calculations for neutral 
Li, Na, and K on helium clusters of size N = 100 — 300 give good qualitative 
agreement with experiment. The doublet structure (^P, ^Pa/o*— ^S'1/2) observed 
in the experimental spectra for Na and K on helium droplets E3 can be seen in the 
PIMC calculations. However, while both experiment and theory show that these 
transitions are shifted to the blue relative to the experimental gas phase values, the 
absolute value of these shifts is in general much more difficult to obtain from theory. 
Due to weak spin-orbit coupling for Li, the doublet splitting is small relative to the 
linewidth, and thus also difficult to resolve. The PIMC absorption spectra for the 
HejvLi system exhibits a weak red shift and a long tail towards the blue, both of 
which are consistent with the experimental spectra. E3 



4.2. Vibrational shifts in infrared spectroscopy of molecules in 
Hcn 

The first spectroscopic experiment made on a doped helium cluster measured the 
infrared absorption spectrum of the octahedral SFg molecule. This low resolution 
spectrum, obtained with a pulsed CO2 laser, revealed that the 1/3 vibrations of SFg 
molecule are red-shifted from the gas phase value by about 1 — 2 cni^^. They also 
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appeared to show that the three-fold degenerate absorption for these vibrational 
modes is split into two peaks. The split peaks were interpreted as implying that the 
SFg molecule resides on the cluster surface where the three-fold degeneracy would 
be expected to be split into parallel and perpendicular modes. However, DMC cal- 
culations made^t that time showed that the molecule should be Jocated at the 
cluster center. £3 This was later confirmed by PIMC calculations Il3 and verified 
by subsequent experimental investigations. These include high resolution spectra 
made with diode lasers EH which showed a single vibrational absorption, red-shifted 
by /S.V = — 1.6 cm~^ and having no splitting of the vibrational degeneracy, and anal- 
ysis of ionization products of SFg-doped helium clusters. Ea Calculations of the spec- 
tral shifts of these triply degenerate intramolecular 1^3 vibrations of SFg were made 
with both PIMC and DMC, ESEj using the instantaneous dipole- induced dipole 
(IDID) mechanism originally proposed by Eichenauer and Le Roy to calculate 
the vibrational spectra of SFg inside argon clusters. PIMC allows the calculation 
of the thermally averaged spectral shift at finite temperatures, while DMC gives 
the ground state, T = K value of the spectral shift. Ref. ^ provides a discus- 
sion of the difficulties in calculating spectral line shapes (and hence extracting line 
widths) from a finite-temperature path integral calculation. The IDID approach is 
taken because the intramolecular vibrational dependence of the He-SFg interaction 
potential is not known, and it is therefore necessary to approximate this. In the 
IDID model, the origin of the spectral shift is assumed to be the dipole-dipole in- 
teraction between the instantaneous dipole moment of the SFg z/3 vibration and 
the induced dipole moments of the surrounding helium atoms. The average shift 
of the 1^3 absorption estimated from the IDID model within PIMC calculations 
are red-shifted, in agreement with experiment, but the magnitudes of the calcu- 
lated shift (Ai/ = —0.84 cm~^) is somewhat smaller than the experimental value 
of Ai/ — —1.6 cm~^. Overall, the agreement of the spectral shift value to within a 
factor of two is quite reasonable, but it is evident that for a proper understanding 
of the spectral shift of SFg inside helium, one needs to also incorporate the contri- 
bution from the repulsive part of the He-SFg interaction that is neglected in the 
IDID model. 

Recently Gianturco and Paesani have calculated vibrationally adiabatic He-OCS 
potentials for various internal vibrational states of OCS, whiciLallow for a more fun- 
damental approach to the calculation of vibrational shifts. l3 These vibrationally 
adiabatic potentials are derived by evaluating the interaction potentials as a func- 
tion of both external coordinates and an internal vibrational coordinate using ah 
initio electronic structure methods, and then averaging over the internal vibrational 
wavefunctions. In particular, they have provided vibrationally adiabatic potentials 
Vbo and Vn that are averaged over the ground state and first excited state of 
the asymmetric stretching motion of the molecule, respectively. The shift of an in- 
tramolecular vibrational mode inside the helium cluster can be estimated within 
an adiabatic separation of the fast intramolecular vibrational mode from the slow 
He-He and He-molecule degrees of freedom. O In this approach, the spectral shift 



Copyright © 2002 by World Scientific, reprinted with permission. 



28 

results from computing the average of the difference between Vbo and Vn over the 
finite-temperature ensemble sampled in the PIMC simulation. Recent PIMC calcu- 
lations for HC39-OCS at T = 0.3 K_hnd a red-shifted asymmetric vibration, with 
the shift of Av — —0.87(1) cm^^. Ed The sign of the shift is in agreement with 
that seen in experimental measurements for OCS made in larger clusters involving 
more than 1000 helium atoms at T = 0.38 K, but its inagnitudc is somewhat larger 
than the experimental value of — —0.557(1) cm^^.lla Detailed analysis indicates 
that these discrepanciffi are likely due to small errors in the vibrationally averaged 
adiabatic potentials. Ed 

4.3. Rotational spectra of molecules embedded in Hcn 

The experimental observation of rotational fine structure for infrared spectra of 
vibrational transitions in the bosonic ^He clusters but not in the corresponding 
fermionic '^He clusters at the operative temperature of T = 0.38 K, EZI led to the 
conclusion that superfluidity is essential for observation of a free rotor-like spectrum. 
This has been explained as a result of the weak coupling of molecular rotations to the 
collective excitations of superfluid He II, compared to the much stronger coupling 
to particle-hole excitations in the Fermi fluid '^He. Ill Consequently, the rotational 
lines are considerably broadened in the fermionic clusters, and the fine structure 
of rovibrational transitions is washed out. This is consistent with the results of di- 
rect calculations of rotational energy levels of ''He clusters containinerotationally 
excited molecules, using zero-temperature DMC-based methods. 13 These 

direct calculations show that the bosonic nature of the ^He is critical in ensuring 
a free rotor-like spectrum of rotational energy levels of the molecule when embed- 
ded in a helium droplet. The corresponding rotational energy levels in fermionic 
helium droplets have not yet been calculated, and would constitute an interesting 
theoretical topic for future study. Path integral calculations have not provided any 
information on the dynamical differences resulting from solvation in fermionic ver- 
sus bosonic helium droplets so far, since fermionic PIMC simulations have not been 
made for these systems. 

A major feature of the rotational spectra of molecules in ^He droplets is the 
appearance of free rotor-like spectra with increased effective moments of inertia. In 
principle, any helium-induced change in the molecular moment of inertia should be 
directly related to a change in the global superfluid fraction, according to Eqs. ( p9| ) 
and (|30| ) (assuming that linear response measures are applicable to the quantized 
rotation of the molecule). Furthermore, a highly anisotropic molecule would be ex- 
pected to result in some anisotropy in the helium response for rotation around 
different axes, yielding anisotropy in the tensor of global superfluid response. tZI 
However, as noted earlier, the global superfluid estimator is relatively insensitive 
to the presence of an impurity and the statistical errors mask small changes. It is 
possible that for significantly larger, and more anisotropic molecules than those the- 
oretically studied to date, e.g., for the planar aromatic molecules such as tetracene 



Copyright © 2002 by World Scientific, reprinted with permission. 



29 

and phthalocyanine that have already been studied experimentally, the global 

superfluid response may be more affected and yield information. For the relatively 
small molecules and complexes studied so far however, it has proven necessary to 
examine the local perturbations of the helium superfluidity in order to develop an 
understanding of the coupling between this and the molecular rotational dynamics. 

The microscopic two-fluid description of the quantum solvation of molecules in 
He that is provided by path integral calculations has led to a detailed analysis of 
the effective moments of inertia of molecules solvated in a bosonic superfluid, and 
hence to a quantitative understanding of the effective rotational constants mea- 
sured in the infrared and microwave spectroscopy experiments. Since the path in- 
tegral calculations carried out to date do not explicitly incorporate the molecular 
rotational degrees of freedom, the connection between the path integral densities 
and the molecular moments of inertia has to be made within a dynamical model. 
Kwon and Whaley have proDOsed a quantum two-fluid model for calculating the 
effective moment of inertia, tllll The main features of this model arc summarized 
below. As will be evident from the assumptions of this quantum two-fluid model for 
superfluid helium response to molecular rotation, it is applicable only to the regime 
of heavier molecules, i.e., those possessing gas phase rotational constants less than 
~ 0.5 cm~^. Excellent agreement with spectroscopic measurements is obtained for 
the two instances in which the He-molecule interaction potential is best known, 
SFg and OCS. i The theoretical values of rotational constant calculated from the 
quantum two-fluid model are 0.033 cm^^ and 0.067 cm^^, for SFg and OCS re- 
spectively, compared with the corresponding experimental values 0.034(1) cm~^ 
and 0.073 cm~^. Draeger et aZJiave recently tested this quantum two-fluid model 
for the linear trimer (HCN)3. EH They also find excellent agreement between the 
predictions of the quantum two-fluid model and the experimentally measured rota- 
tional constant. While the HCN monomer lies in the regime of light molecules, the 
trimer is sufficiently massive to fall within_the heavy regime, possessing a gas phase 
rotational constant of Bq — 0.015 cm~^. E3 

The quantum two-fluid model of Kwon and Whaley is a microscopic two-fluid 
continuum theory for the spectroscopic response of a molecule rotating in super- 
fluid He. It is to be distipguished from the phenomenological two-fluid theory of 
Landau for bulk He II. Ej The Kwon and Whaley model examines the helium re- 
sponse to rotation of an embedded molecule that, starting from the local two- fluid 
decomposition of the molecular solvation density into non-superfluid and superfluid 
components. Unlike the Landau two-fluid theory for bulk He II, it does not make 
a two-fluid decomposition of the current densities, but deals only with decompo- 
sition of the helium density near an impurity on an atomic length scale, into a 
non-superfluid component induced by the molecular interaction and the remaining 
superfluid. This constitutes a signiflcant difference between the well-known phe- 
nomenological theory for bulk, homogeneous He II, and the microscopic two-fluid 
model for molecular rotational dynamics in an inhomogeneous superfluid solvation 
situation. For the remainder of this section we shall interchangeably use the terms 
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two-fluid model, microscopic two-fluid model, and quantum two-fluid model to refer 
to the Kwon/Whaley model. 

The starting point for the quantum two- fluid model for helium response to molec- 
ular rotation is the local two-fluid density decomposition of the molecular solvation 



density that results from path integral calculations. As described in Sec. 3.4, con- 
sistent evidence for the existence of the local non-superfluid density in the first 
solvation shell around the molecule, induced by a strong molecular interaction with 
helium, has now been obtained from two different estimators of the local superfluid 
response. The second feature of the model is the assumption of adiabatic following 
of some or all of the solvating helium density with the molecular rotation. Adiabatic 
following means that the helium density follows the molecular rotational motion. 
Quantitatively, complete adiabatic following of the helium density would imply that 
when viewed in the rotating molecular frame, the helium density appears station- 
ary. Thus, in the molecular frame it is independent of rotational state. This applies 
to both classical and quantum descriptions of the molecular rotation. In a classical 
description the helium density is analyzed as a function of continuous molecular 
rotation frequency, while for a quantum description it is analyzed as a function of 
quantum rotational state of the molecule. 

The accuracy of the adiabatic following assumption, as well as quantification 
of the extent of adiabatic following_bi!_hclium, has been the subject of several 
studies by Whaley and co-workers. U'ESe^I Within a classical description of molecu- 
lar rotation, Kwon et al. provided a criterion for adiabatic following, namely that 
the kinetic energy of rotation associated with a particular helium density (total, 
non-superfluid, or supfirfluid) be less than the potential energy barrier to rotation 
around the molecule. □ This criterion is applicable to densities deriving from any 
number of helium atoms, and allows simple estimates using either barriers to rigid 
rotation, or barriers to adiabatic motion between potential minima associated with 
different molecular orientations. Application of this criterion to the molecules OCS, 
SFg and HCN, for which the molecule-helium interaction potentials are very well 
characterized, showed that for both OCS and SFg it is energetically feasible for the 
entire helium density to adiabatically follow the molecular rotation. However, for 
the lighter HCN molecule, it is not energetically feasible for any density component, 
whether non-superfluid, superfluid, or total, to adiabatically follow the molecular 
rotation. In a classical sense, the molecular rotation is then too fast for the helium 
to follow. The consequence of this lack of adiabatic following is that the helium 
density distribution is more diffuse when viewed in the molecular frame. This can- 
not be seen directly in the PIMC densities, since molecular rotation is not included 
in these. However, it can be seen directly in diffusion Monte Carlo calculations of 
excited rotational states of the molecules in He clusters. In these calculations jnade 
with an importance sampling algorithm for rotational degrees of freedom, cil the 
helium density (or wave function) is projected into the rotating molecular frame 
and compared with the corresponding density (wave function) from a calculation 



performed without molecular rotation. 



Explicit analysis of the dependence on 
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rotational state can also be performed, although comparison between rotating and 
non-rotating cases is already very revealing. The original application of this anal- 
ysis showed that the extent of adiabatic following decreases for lighter molecules, 
with the helium density in the molecular frame becoming more diffuse as the rota- 



comparison may be quantified by evaluation of a quality factor Q that measures 
changes in the ratio of densities along directions corresponding to strong and weak 
binding, as a function of molecular rotational state. B Complete adiabatic following 
is measured by Q = 1, provided the molecular interaction potential is anisotropic. 
(For an isotropic interaction with helium, adiabatic following is not applicable, and 
Q = 1 by definition.) Application to the series of molecules, OCS^Fg, and HCN, 
shows that Q - for HCN, and Q ^ 0.7 for both OCS and SFe. This confirms 
the prediction of the energetic criterion for HCN, i.e., there is negligible adiabatic 
following around this molecule. The Q-value results for the heavier molecules are 
quite significant, implying that the extent of adiabatic following is not complete, 
even for the most strongly bound case of a single He atom attached to SFg. i So 
only a fraction of the helium density can adiabatically follow the molecular rotation, 
even for a heavy, strongly bound molecule. 

The next stage of the quantum two-fluid model is to consider the consequences 
of adiabatic following for both the local non-superfluid and local superfluid density 
around a dopant molecule. These two density components show very different re- 
sponse to adiabatic following, deriving essentially from the different spatial extent 
that results from their corresponding underlying exchange permutation paths. The 
molecule-induced non-superfluid density is localized close to the molecule, within 
the first solvation layer, and is composed of very short permutation exchange paths. 
In order to satisfy adiabatic following, such a localized density must rotate rigidly 
with the molecule. There is no other obvious way in which a density that is spatially 
localized within a few angstrom can remain constant in a rotating molecular frame. 
This results in an increment of moment of inertia from the local, molecule-induced 
non-superfluid that is given by 



For the heavy molecules SFe, OCS, and the linear trimer (HCN)3, the PIMC values 
for A/„ amount to 100%, 90%, and ~ 81% of the corresponding experimentally 
observed moment of inertia increments. A/. It is interesting that for the highly 
symmetric SFe molecule, a very similar result (AJ„ ~ 98% of A/) is obtained from 
calculations with only an isotropic interaction potential. While there is no adiabatic 
following with an isotropic interaction and hence no mechanism for rigid coupling 
of the non-superfluid helium density to the molecular rotation, the high symmetry 
of the octahedral SFg molecule nevertheless results in the integrated non-superfluid 
density in the first shell being very similar in anisotropic and isotropic calcula- 
tions. In fact, the finding that the isotropic non-supcrfiuid density could account 



tional constant of the molecule increases. 



Kwon et al. showed recently how this 




(4) 
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quantitatively for AL was obtained prior to calculations of the anisotropic local non- 
superfluid density. While this result did not have the theoretical justification of 
rigid coupling as a result of adiabatic following at that time, it was the first indica- 
tion that a local two-fluid description was dynamically relevant and prompted the 
application of an microscopic AncLmnikashvili analysis of experimental rotational 
spectra for the case of OCS in He. 113 

In contrast to the local non-superfluid density, the superfiuid density, while also 
modulated around the molecule, is not restricted on an angstrom length scale within 
the quantum solvation structure. By its very definition, consisting of long exchange 
paths, the superfiuid density extends far from the molecule. Thus the equation of 
continuity can applied to this density over long distances. Kwon et al. have shown 
that for a classical molecular rotation, determined by a continuous frequency w, the 
condition of adiabatic following, if satisfied, can be combined with the equation of 
continuity to eliminate the explicit time dependence of the density and to arrive at 
an equation for the superfiuid velocity: 13 

V • [p,(r,i)v,(r,t)] = Vp,(r,t) -{ujxv). (5) 

The irrotational nature of a superfiuid may be used to replace Vs by (?i/m4)VM(r, t), 
to arrive at a second-order partial differential equation for the superfiuid velocity 
potential u{t): 

V • [ps{v)Vu{v)] = (^) Vp,(r) . {u X r). (6) 

This equation, discussed in detail in Ref. ^, was first proposed in before full 
anisotropic superfiuid densities in three dimensions were calculated. c3 

A similar equation was recently presented by Callegari and co-workers, to- 
gether with the somewhat different assumption that the entire local solvation den- 
sity is superfiuid. Solution of these hydrodynamic equations leads to a hydrody- 
namic moment of inertia increment A/^ that is derived from the excess fluid kinetic 
energy associated with the flow pattern of v^. Callegari et al. solved the hydro- 
dynamic equations for several linear or rod-like molecules for which the equations 
becomelwo-dimensional, using total densities derived from density functional calcu- 
lations. The full solution for a molecule showing true three-dimensional anisotropy 
was made recently for SFg using PIMC densities. El Draeger et al. have applied the 
hydrodynamic treatment to the (HCN)3 trimer, using their PIMC densities and 
also assuming the total density to be superfluid, for the purpose of comparison. 
It appears that very different results are obtained for different molecules within 
the hydrodynamic treatment. For octahedral SFg, the value of A//i is small, irre- 
spective of whether the total density or superfiuid density is used as input to the 
hydrodynamic calculations (6% and 9% of A/, respectively). For the linear (HCN)3 
trimer, Draeger et al. find an upper bound of AJ/j ~ 0.7A/, when the total density 
is assumed superfiuid (A/;; = 850 amu A^, compared with an experimental value l3 
of A/ = 1240 amu A^). The calculations of Callegari et al. for rod- like molecules 
yielded between A//; ~ 67% and 98% of the experimentally measured increments 
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A/. These studies differed from those for SFg and (HCN)3 in that input densities 
were obtained from density functional calculations rather than from PIMC, in some 
cases using simple estimates from pairing rules to construct interaction potentials 
when no empirical or ab inito potentials were available. 

The hydrodynamic treatment of the local superfluid density derived from PIMC 
has a number of questionable aspects, i Firstly, the treatment of the molecular 
rotation as a classical rotation characterized by a continuous frequency w must be 
reconciled with the intrinsic quantized nature of spectroscopic transitions between 
quantum rotational states. The response to classical rotation necessarily gives rise 
to angular momentum generatiQii_in the superfluid, analogous to the rotation of 
bulk superfluid in a superleak. EilEa Kwon et al. have calculated the angular mo- 
mentum generation by absorption of a photon within a semiclassical analysis, and 
shown that significant values of Alh result in large fractions of the photon angular 
momentum being transferred to the superfluid density component. This contradicts 
conclusions of a number of zero-temperature DMC-based calculations that indi- 
cate there is nngl igible transfer of angular momentum to the fluid on rotational 
excitation. ElcilO Kwon et al. resolved this by adding quantum constraints to the 
hydrodynamic formulation, and concluding that violation of these indicates invalid- 
ity of the hydrodynamic contribution. This in turn may derive from lack of complete 
adiabatic following, for which considerable evidence now exists, as outlined above, 
or from the intrinsic lack of applicability of hydrodynamics to the motions of a 
superfluid on the atomic length scale. An indicator of this breakdown is the fact 
that the solutions to the hydrodynamic equations with density inputs of atomically 
modulated helium solvation densities around an embedded molecule, show varia- 
tions over length scales of 1 to 2 A. i-El Such variations on a distance comparable to 
or less than the coherence length ^ of helium imply that a hydrodynamic solution is 
at its limits of validity here, at best, and should be interpreted with great caution. 

The overall conclusions of the quantum two-fluid model for the response of he- 
lium to rotation of an embedded molecule are thus that the primary contribution 
to the increased molecular moment of inertia is a rigid coupling to the local non- 
superfluid density in the flrst solvation shell. This yields 100%, 90%, and ~ 81% of 
Al for the heavy molecules SFg, OCS, and (HCN)3, respectively. The accuracy of 
these estimates is dependent on the accuracy of the underlying molecule-helium in- 
teraction potentials. In contrast, there appears to be negligible contribution from the 
superfluid, whose response must be restricted by angular momentum constraints. 
This is consistent with the flndings of only partial adiabatic following of the total 
helium density. It appears reasonable that only the non-superfluid density adiabat- 
ically follows the molecular rotation, while the superfluid density, which is defined 
over much longer length scales, cannot effectively adiabatically follow. For heavy 
molecules, this two-fluid model provides a complete dynamical picture. 

For light molecules such as HCN, the zero-temperature calculations have shown 
that adiabatic following is questionable even for the non-superfluid density. Conse- 
quently, in this situation the two-fluid model cannot be used to estimate effective 
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moments of inertia. At this time, the zero-temperature DMC-based direct calcula- 
tions of rotational energy levels of doped clusters provide the only route to micro- 
scopic theoretical understanding of spectroscopic rneasurements of rotational tran- 
sitions for such light molecules in helium droplets. C3 This will hopefully change in 
the future, when molecular rotational motions are explicitly incorporated into the 
PIMC. 



5. Conclusions and future directions 

The path integral approach has provided a powerful theoretical tool for investigat- 
ing the superfluid properties of finite helium droplets. Path integral Monte Carlo 
calculations have shown that these systems constitute nanoscale superfluids and 
offer a unique route to probing the structure and response of a Bose superfluid 
on a microscopic length scale. They also provide examples of inhomogeneous super- 
fluid density, with the unique feature that this inhomogeneous, nanoscale superfluid 
density can be probed by molecular and atomic dopants. 

The microscopic calculations show that such impurities introduce a local quan- 
tum solvation structure into the otherwise smoothly varying helium density. The 
numerical path integral Monte Carlo method has allowed this quantum solvation 
structure in a superfluid to be analyzed in terms of the boson permutation exchange 
properties, and conversely, the effect of the molecular interaction on the superfluid 
to be quantified. PIMC calculations show that a strongly bound impurity induces a 
non-supcrfluid density in the first solvation shell, whose extent is determined by the 
strength of the molecular interaction. Similar conclusions are derived from analysis 
of the permutation exchange path lengths into short, strongly localized paths, and 
long, delocalized paths, and from decomposition of the linear response estimator 
for global superfluidity. The existence of this local non-superfluid component in the 
solvation layer around microscopic impurities therefore seems to be a general fea- 
ture of molecular solvation in superfluid He clusters. Response of this local two- fluid 
density to the rotation of a molecular impurity gives rise to increments in the molec- 
ular moment of inertia, but does not otherwise modify the effective free rotation of 
the molecule in the superfluid. 

These path integral studies of doped helium droplets open the way to study 
of several intriguing questions. One is the effect of the molecular rotation on the 
quantum solvation structure in the superfluid local environment. As noted in this 
article, all PIMC studies to date have not explicitly incorporated the molecular 
rotational degrees of freedom. We now know from zero-temperature calculations 
of the quantum rotational excitations that the molecular rotation does result iiL_a 
smearing out of the angular anisotropy in the quantum solvation structure. Ij'OLiI 
This implies less than perfect adiabatic following of the helium density, even in the 
rotational ground state. Given the significance of the adiabatic following assumption 
for models of the helium response and hence for the analysis of rotational spectra 
of doped molecules, developing a direct route to the solvation density around a 
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rotating molecule is highly desirable. This can be done by incorporating the molec- 
ular rotational kinetic energy in the path integral representation. A key question 
with spectroscopic implications is then how the local two-fluid density decompo- 
sition is modified. We noted earlier that the moment of inertia increment of the 
non-supcrfluid density around SFg is approximately independent of the anisotropy 
of the interaction potential. This suggests that even if the two-fluid densities are 
modified with rotation, becoming less anisotropic, the effective moment of inertia 
of the molecule in ^He will be unchanged. This remains to be verified. 

A second direction departing from the analysis of molecular solvation structure 
in a superfluid is the investigation of localization of helium atoms and their removal 
from the superfluid state, as a function of the binding to organic molecules of in- 
creasing size. In the study of benzene, the key feature responsible for the localization 
phenomenon was identified as the strong and highly anisotropic interaction of he- 
lium with the TT-electron system. Systematically varying the extent of the 7r-electron 
system by going to larger planar, polyaromatic molecules will allow the transition 
from a nanosubstrate to a microsubstratc that begins to mimic a bulk solid surface 
to be investigated. We expect that the "inert" layers familiar from studies of thin 
films of helium on graphite will evolve from these localized atoms, but the manner in 
which this happens will depend on the role of lateral confinement and permutation 
exchanges in the presence of an extended 7r-electron system. 

A third, novel direction is provided by extension of these ideas to nanoscale 
clusters of molecular hydrogen, H2. In its rotational ground state, the II2 molecule 
is a boson, but bulk superfluidity is preempted by the occurrence of the triple point 
at T = 13.6 K. However, finite-size and reduced dimensionality systems are offer 
ways of bypassing this solidification of hydrogen by allowing lower densities and 
thereby moderating the effects of strong interactions. Path integral calculations have 
already been used in several instances in the search for a superfluid state of molecular 
hydrogen. Thus, very small finite clusters of (H2)Ar (N < 18) have been shown with 
PIMC to be not only liquid-like but also to show a limited extent of superfluidity. E3 
Two-dimensional films of hydrogen have been shown to allow a stable superfluid 
phase at low temperatures provided that an array of alkali atoms is co-adsorbed, 
providing stabilization of a low density liquid phase. Ell Given these low-dimensional 
antecedents, it appears possible that a relatively small solvating layer of hydrogen 
wrapped around a molecule might also show some superfluid behavior. Path integral 
calculations are now in progress to examine the extent of permutation exchanges in 
cycles around different axes of a linear molecule wrapped with variable numbers of 
H2 molecules. Such studies will provide microscopic theoretical insight into the 
quantum dynamics underlying recent spectroscopic experiments showing anomalies 
in the molecular moment of inertia that are consistent with a partial superfluid 
response of the solvating hydrogen layer. E3 

In summary, the path integral Monte Carlo approach provides a unique tool for 
analysis of these degenerate quantum systems in finite geometries and with chem- 
ically complex impurity dopants. The insights into nanoscale superfluid properties 
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that have resulted, and the interplay between physical and chemical effects afforded 
by calculations on doped helium droplets offer promise of new opportunities for 
analysis and manipulation of supcrfluid at the microscopic level. 
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